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Abstract 

Ant algorithms are inspired in real ants and the main idea is to create virtual ants that travel 
into the space of possible solution depositing virtual pheromone proportional to how good a specific 
solution is. This creates a autocatalytic (positive feedback) process that can be used to generate 
automatic solutions to very difficult problems. In the present work we show that these algorithms 
can be used coupled to tight-binding hamiltonians to design conducting polymers with pre-specified 
properties. The methodology is completely general and can be used for a large number of opti- 
mization problems in materials science. 

PACS numbers: 
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INTRODUCTION 



The design of materials with pre-specified properties represents a complex and very diffi- 
cult problem. Recently, many groups have attempt to use artificial intelligence or automatic 
methods to design new materials [l, 2], ^,0, 

Real ants consist in social insects that, while individually they can not be considered 
very intelligent creatures, but as a colony they display remarkable capability of solving very 
complex optimization problems. This kind of collective intelligence of the swarm {(], 0] can be 
considered an emmergent property of the system since it appears to be due to the high level 
of interaction among the agents (ants). The problem of finding food in complex environments 
is similar to some hard computational problems that can not be handled by conventional 
searching methodologies. Examples of these problems are present in distinct areas such as 
protein folding, logistics and the design of materials with predetermined properties. 

Based on the capabilities of real ants to solve hard optimization problems (such as finding 
best or quasi-best paths connecting a food source and the ant nest) a new class of meta- 
heuristic algorithms, named ant system (ANT), was recently developed B, H, fl 3- The 
mechanism that allows the ants to establish optimized paths is based on the fact that when 
ants move they deposit chemical markers (pheromones) , which defines specific trails to be 
travelled. The ants tend to move along trails where the pheromone concentration is higher. 
This kind of process is known as an autocatalytic one due to the positivie feedback. For 
instance, if two ants establish two different routes connecting food sources and the nest, 
the process of going back and forth along these different routes will increase the amount 
of pheromone in the most best one (Fig. la-lc), since it will be more walked in the same 
amount of time. Higher concentration of pheromone in one specific trail increases the proba- 
bility of another ant moving along this trail and so on, generating the autocatalytic process. 
This is equivalent to a computational parallel processing probing the space of good solutions 
in a fast and efficient way. 

These aspects of the ant organization were computationally translated into the ant al- 
gorithm. The idea is to create artificial (virtual) ants that 'travel' in trails of the phase 
space of possible solutions. These ants can 'communicate' through the virtual 'pheromones' 
deposited on these trails, with the deposition rate proportional to the effectiveness of the 
tentative solution (Fig. lb). As in the case of real ants the positive feedback in the vir- 



2 






(1) 


!A.B+dB] 


d 




- — i 


► — » 






■ 




J I2] 




13] 


|A,9-dB) 





FIG. 1: Schematic view of the ant algorithms. See text for discussions. 



tual space of solutions is capable of generating good solutions that 'autocatalitically' evolve 
in time, in general being able to solve very hard optimization problems 

fl flS,B|. The 

ant algorithm implicitly possessing abilities of 'memory' and 'learning' can outperform even 
recent bioinspired class of genetic algorithms [lol ]. 

Genetic algorithms consist in a class of bioinspired algorithms that were originated from 



the studies of J. Holland in the 1970s 



10] . The metaphor underlying this method is that of 
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natural evolution, whose concepts are followed in a very simple scheme that allows computers 
to evolve automatic solutions over time. The method is constructed based on the evolution 
of populations using processes analogues to the biological ones, namely the mutation and 
the crossover (variation, selection, and inheritance). Genetic algorithms are very effective 
to obtain good solutions from "scratch" , but they are less efficient to improve on good ones. 
The reasons are the same as in real biology. As the available genetic space starts to exhibit 
little variations, the mechanisms of evolution (in special selection) are no longer very effective 
and new solutions must heavily depend on the mutations, which in general are effective only 
over long periods of time. Although there are some computational tricks to overcome these 
difficulties and to improve on the algorithm performance, sometimes they are very difficult 
to implement numerically. Thus, it would be desirable to have alternative methodologies. 
The class of ant algorithms can be one of these alternatives. 

In the present work we adapted ant algorithms to solve the problem of designing con- 
ducting polymers. Polymers constitute a very large and important class of materials that 
exhibit a rich variety of mechanical and electronic properties, varying from soft to very hard, 
from insulators to high conducting materials [111 ]. Due to the richness of the carbon atom 
reactivity an almost infinite number of new structures is possible. This makes the system- 
atic search for new structures impossible and the trial and error approach has been the rule. 
In this sense the use of combinatorial and/or artificial intelligence methods to assist in the 
design of new or to improve on the existing materials would be very usefull. In order to show 
that ant based algorithms can be an effective tool to do this we illustrate its application to 
the complex problem of designing polymeric alloys with pre-defined electronic properties. 
We considered the case of the polyaniline (PANi) family, one of the most studied conduct- 
ing polymers l^ . PANi encompass a family of compounds where nitrogen atoms connect 
six-membered carbon rings of benzenoid or quinoid character. They can have three differ- 
ent monomeric units (PAN1,PAN2,PAN3, see Fig. 2) depending on its degree of oxidation 
and/or protonic doping. The relative concentration of these units in the polymeric chain 



define its conductivity, varying from insulator to highly conducting material [12J, I3I Il4l. Il5|. 
Considering that, in general, each polymeric chain contains hundreds of PANI, PAN2, PAN3 
units [12I a huge number of possible configurations varying the relative concentration of the 
different units is possible. What are the concentrations of PANI, PAN2, PAN3 to produce 
the most conductive alloy (metallic states) ? The metallic state is characterized by a zero 
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gap value and for delocalized states at the Fermi level. We coupled our based ant algorithm 
to the simple Huckel hamiltonian (one of simplest tight-binding models) 16| to solve this 
problem. We were able to automatically generate very good solutions. The methodology is 
completely general and can be coupled to sophisticated hamiltonians (even ab initio ones) 
and can be applied to design other classes of materials such as metals or oxides. 
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FIG. 2: Structure of the three monomeric units used in the calculations: PAN1, PAN2, PAN3. 



METHODOLOGY 

Ant System 

Our adapated ant algorithm is based on two principles: positive and negative feedback, 
which are associated with the creation and destruction of trails, respectively. The trail cre- 
ation occurs when an ant walks between two points on a potential surface P(x, y), depositing 
an amount of pheromone that is proportional to the quality of the solution at the end point. 
This process is described by the trail function: 

T(x 1 ,yi 1 x 2 , y 2 , t) = T{x u yi, x 2 , y 2 , t - 1) + P(x 2 , y 2 ) (1) 

where (x\, y\) and (x 2 , y 2 ) are the the starting and end points (see Fig. Id), respectively, and 
t is the current time step. The potential function P(x, y) is the mathematical translation 
between this method and the problem to be solved. This function is the responsible for the 
positive feedback of good solutions since better points will generate better (more walked) 
trails. The trail is the most important interaction feature in this method since it is the only 
communication mechanism among the ants. 
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The contraposition to this pheromone accumulation in certain points is given by the 
regulatory negative feedback contribution, introduced as a constant evaporation rate term, 
p < 1, acting over the trail function: 

T(x lt yi, 22, 2/2, t) = pT(xi, yi,X2, y 2 , t). (2) 

This term is the responsible for the removing of the trails that are not being walked, 
since in these trails there is no pheromone redeposition. This is a strong selection acting 
over trails leading to bad solutions, contributing to the restriction of the circulation in the 
surface and to the saving of time in the searchings. 

This pheromone deposition/evaporation process has the only objective of guiding the ant 
movement choices when faced with the options given by potential surface. The choices are 
quantified by an evaluation function: 

E{xi,y 1 ,x 2 , 2/2, t) = aP(x 2 ,y 2 ) + bT(x 1 ,y 1 ,x 2 ,y 2 ,t) + cR(t) (3) 

which is calculated for each of the possible neighboring final points connected to the current 
point occupied by the ant. The tuning coefficients a, b and c are used to calibrate the relative 
weights between the contributions. R(t) is a random term which represents a kind of noise 
or uncertainty in the information acquired from the environment, leading to the appearance 
of small fluctuations which are relevant in situations where two different evaluation functions 
have very close values. The point with the highest value for the evaluation function is the 
one which will be the next to be visited. 



Polymer Problem 



As mentioned before our objective is the determination of the best set of relative concen- 
trations (x(PANl),y(PAN2),z(PAN3)) (see Fig. 2) that produces the polymer with better 
electrical conductivity. One of the objectives here is to demonstrate that ant algorithms can 
solve this kind of problem in a fast and efficient way without extensively probing the space 
of possible solutions. Thus, we started analyzing a size problem that can be solved by sys- 
tematic search and then we compare the 'exact' solutions with the ones obtained using the 
ant algorithms. In this sense we considered PANi chains containing 100 (equal or different) 

ectronic structure we used the simple Hiickel (tight-binding 



units. In order to obtain their e 
model) hamiltonian I3I, 14, Ij], 



16J: 



h = 5>|i>(*i+ E /woi, ( 4 ) 

where a, and represents the Coulomb and the hopping integrals, respectively [l6(. The 
wavefunctions are written as a linear combination of n atomic 7r-orbitals 

n 

I*) =S C «l0i)> ( 5 ) 

The degree of electronic wavefunction delocalization can be estimated from the inverse 
participation number (IPN) [itJ given by: 



IPN = ^t 1 ' • (6) 

where q are the coefficients of the highest occupied molecular orbital(HOMO) (Fermi 
level) . 

IPN can assume values from (maximum delocalization, the wavefunction is spread over 
the whole polymeric chain) up to 1 (maximum localization, the wavefunction is localized 
over a single atom). From the IPN and gap (energy difference between the frontier (last 
occupied - first unoccupied) molecular orbitals) we can characterize the metallic states. In 
order to be a good conductor the gap value should be smaller than the room temperature 
thermal energy (25 meV) and the IPN value close to zero. We need do determine the x,y, 
and z values that satisfy the gap and IPN conditions. The problem can be simplified if we 
use the constraint relation x + y + z = 100, which reduces the space of possible solutions from 
a tridimensional (x, y, z) to a two-dimensional (x,y surface) one. If the polymeric chains are 
long enough the relative concentration of the different units is more important than their 

nnn 

specific position in the chains, considering that they are randomly formed [13|, H4], 115[ • In 
this case we would have a grid of 100x100 of possible configurations that by the symmetry 
imposed by our contraint relation is then reduced to a grid of 5151 points. Each point 
{xkiUk) (k = 1, • • • , 5151) in the grid represents a possible solution (with specific gap and 
IPN values) of our problem associated with a polymeric chain containing 100 monomeric 
units distributed accordingly to the relative concentrations (xk,Vk^k)- We then proceed to 
calculate the gap and IPN values for each one of these 5151 configurations. The calculated 
space of solutions can be better visualized through a potential function of the grid points 
P{x k ,y k ): 



P{x <« »> = log (gap^/xipn^)) (7) 

We used a /og function due to the large variations found for neighboring points. The 
surface topography of the grid points calculated with this function is presented in Figure 
[3j The observed triangular shape is a consequence of the used constraint relation. Higher 
peaks are indication of metallic signatures. The visual complexity of this profile gives an 
idea of the difficulty of the problem for standard optimization algorithms that in general fail 
for problems where many local minima and multiple local maxima are present. 




% PftH 

FIG. 3: Top: Contour plot of the potential grid for PANi chains containing 100 units as a function of 
the percentage of units x(PANl) and y(PAN2); Bottom: Zoomed contour plot for the initial region 
of the grid. The corrugation of the grid and the existence of many local maxima are indications of 
the complexity of the potential profile. 
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Adapted Ant Algorithm 



Our adapted ant algorithm can be coupled to the Huckel hamiltonian [l6| to carry out 
the automatic search for the metallic structures. The general algorithm fluxogram has the 
following form: 



To generate the starting positions 
do i=l, number of cycles 

do j=l, number of steps 

do k=l , number of ants 
do 1=1,4 (neighbours) 

if the next point is permitted 

if the electronic properties of this point were 
not calculated 

Calculate the electronic properties 
Calculate the potential function 
Calculate the evaluation function 
end do 
Choose the next point 
Move the ant 

Update the pheromone trail 

Check if the stop condition was achieved 

end do 

Multiply the trail function by the evaporation term 
end do 

To generate the ants positions for the next cycle 
end do 



The algorithm is started with the random sorting of the ant starting positions. There 
are two temporal variables: the cycles and the steps. The cycle is a fixed number of steps 
after what the ants are sorted again over the grid. Given a point in this grid the ant has 
four possible directions (up, down, left, right) (see Fig. Id). In order to avoid closed cycles, 
each ant stores the sequence of points visited and these are prohibited to be visited again by 
the same ant. There are other two prohibitions for the ant movements: if there is another 
ant in the final point (since two ants are not allowed to have the same coordinates at the 
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same time) and if the final point is outside the region defined by the constraint relation. 
Everytime when the electronic parameters of a given point are calculated the information 
is stored in a databank accessible to all ants. This prevents redundant calculations and 
optimize the searching time. If no improvement on the best solution found is obtained after 
the maximum number of allowed cycles the simulation is ended. 

The ant system was tested in 100 runs with different initial positions for the ants and 
different seeds for random number generator. We have used 50 ants, 100 steps, and 5 cycles. 
The tuning parameters were a = 1.0, b = 10, c = 0.1, and p = 0.95 (see Eqs. 2 and 3), and 
were chosen based on previous testing calculations. The number of ants is approximately 
1% of the total number of grid points, condition that we established as a rule of thumb. The 
algorithm seems to be very robust. In test simulations the obtained results were not very 
sensitive to the tuning parameters. 

The strategy to quantify the efficiency of the method is not based on the number of 
steps but on the average number of electronic calculations iV over the independent 100 runs 
necessary for finding a specific grid point, since this is the main time consuming operation in 
the simulations. As mentioned earlier more then one ant can pass each point, but this does 
not imply in a repeated electronic calculation, since that for the grid point probed the gap 
and IPN values were already stored. N in combination with the percentage of grid points 
probed in the set of simulations (Nf), can provide a good measure of the efficiency of the 
methodology. 

RESULTS AND DISCUSSIONS 

In Table [I] we present a summary of the obtained results. The composition (number of 
PAN1, PAN2, PAN3 units), gap and IPN values for the 10 best solutions obtained from 
the systematic search are displayed. All these points have effectively a zero gap value since 
these energies are smaller than the thermal energy at 300 K (25 meV). The associated IPN 
values are very small, thus indicating a significant electronic derealization at the Fermi 
level. Thus, all these 10 points satisfy the required conditions for metallic states. As the 
gap values are pratically zero and the IPN difference between the best and the worst case 
is only 5.8 x 10~ 4 eV, from a practical point of view these ten best grid points are basically 
equivalent. 
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TABLE I: Values for GAP, IPN, P, N f , and N for the ten best grid points. 



(x k ,y k ,z k ) GAP(10- 6 


eV) IPN(10- 2 


eV) P(x k ,y k ,z k ) 


Nf (%) 


N 


(93,5,2) 


0.67 


1.072 


18.745 


88 


1768 


(91,6,3) 


0.74 


2.009 


18.025 


79 


1557 


(95,2,3) 


19.58 


0.167 


17.236 


88 


1547 


(9,90,1) 


7.24 


0.618 


16.922 


38 


999 


(14,84,2) 


23.46 


0.259 


16.622 


82 


1431 


(40,59,1) 


4.97 


1.302 


16.553 


47 


1148 


(41,58,1) 


4.97 


1.302 


16.553 


58 


1265 


(38,61,1) 


4.97 


1.302 


16.552 


79 


1433 


(39,60,1) 


4.97 


1.302 


16.552 


79 


1489 


(42,57,1) 


4.99 


1.130 


16.548 


63 


1350 



The results obtained from the ant algorithm: the potential function P(x k ,y k , z k ), the 
percentage of the 100 runs where the specific grid point was located (Nf) and the average 
total number of points (N) probed are also displayed in Table [I] 

We can see from the Table that the majority of the 10 best points were located in 
the search runs and in average less thant 30% of the space of solutions was probed. The 
best 3 solutions were located in about 80% of the runs. This is a remarkable performance 
considering the existence of multiple minima and maxima and that these points are basically 
equivalents (very close gap and IPN values). 

An example of the trail evolution of one simulations is presented in Figure 4. The frame 
sequence starts with the initial random ant distribution (Figure 4a). In the first steps the 
pheromone term contribution (see Eqs. 2 and 3) is not yet important since there is pratically 
no trail overpositions. The situation changes when the ants start to find 'weak' trails formed 
by single ants. Due to the importante of pheromone contribution to the evalution function 
(b=0, see Eq. 7), there is a tendency of the of ant to follow these starting trails (Figure 4b). 
In fact as more ants find these trails, larger will be the amount of pheromone deposited and 
greater the importance of these paths (Figure 4c). The global evolution of the system goes 
in a direction where the interconnected trails form a network of preferential paths accessible 
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to all ants, and the positive/negative feedback established (Figure 4d). 
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FIG. 4: Evolution of the trail distribution with time. (A) Initial random sorting; (B) Evolution 
before the start of trail overlap (20 steps); (C) Formation of the network due to the trail overlap 
(90 steps); (D) Final set of trails formed (25 steps of cycle 4). 
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However, this is not a static network since the pheromone evaporation continuously acts 
in the removal of trails. The continued existence of one trail in the long range depends on 
the quality of the points contained in it. This effect is observed when the final distribution 
of pheromone is plotted superimposed with some of the 10 best points shown on Table [B 
Figure shows this final distribution at the end of the calculation shown in Figure HI 



FIG. 5: Pheromone distribution map at the end of the calculation of Fig. |H The five marked 
points were selected from the top 10 list (Table H|). 

It can be noticed that all points showed are in trails that survived after evaporation. This 
demonstrates that the evaporation (negative feedback) acts mainly in removing bad points 
where the ant rarely travels. This also shows why the method is so effective: the trails 
passing through some of the best points are continuously reinforced (positive feedback) and 
the random term (last term of Eq. 3) warrants the possibility of finding yet better solutions, 
as the case of real ants. A better idea of the dynamics of these processes can be obtainded 
from the movie showing an animated sequence of the virtual ant walkings into the space of 
solutions 

SUMMARY AND CONCLUSIONS 

The design of materials with pre-specified properties represents a complex and very dif- 
ficult problem. Ant algorithms are a new class of bioinspired algorithms capable of solving 
very difficult optimization problems. We present here a new methodology for automatic de- 
sign of polymeric alloys with pre-especified properties, in our case we are looking for metallic 
signatures. We coupled an adapted ant algoritms to the simple Huckel hamiltonian (one of 
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simplest tight-binding methods) in order to determine the composition of polymeric alloys 
based on polyaniline polymers capable of exhibiting metallic behavior. Our results show 
that the ant algorithm is efficient and well suited for the optimization problem involving 
the design of polymeric alloys. We were able to obtain very good solutions for the PANi 
alloys, probing less than 30% of the space of solutions. Preliminary results using the same 
methodology for ternary and quartenary alloys of very long polymeric chains, where the sys- 
tematic search is computationally cost prohibitive, are very encouraging. Our methodology 
is hamiltonian independent, it could be used with other hamiltonians, even ah initio ones. 
It could be applied to other classes of materials, such as metals or oxides. We hope the 
present work to estimulate further works along these lines. 

The authors acknowledge financial support from IMMP/MCT, THEO-NANO and the 
Brazilian agencies FAPESP, CNPq and CAPES. 
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